{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a218fde8",
   "metadata": {},
   "source": [
    "# Optimization Gone Wrong and How to Fix it!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd825a91",
   "metadata": {},
   "outputs": [],
   "source": [
    "%pip -q install gurobipy==13.0.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a642764",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import gurobipy as gp\n",
    "from gurobipy import GRB"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc6fdbb6",
   "metadata": {},
   "source": [
    "## Getting to Know a Parameter with Gurobot\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6958426b",
   "metadata": {},
   "source": [
    "Here is a variation of the `quadratic assignment problem (QAP)`. Below is a formulation. We really won't need it but it's here for fun."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8beed826",
   "metadata": {},
   "source": [
    "#### QAP Mathematical Formulation\n",
    "\n",
    "##### Sets and Indices\n",
    "- $i, k \\in \\{0, 1, \\ldots, n-1\\}$: Facilities\n",
    "- $j, l \\in \\{0, 1, \\ldots, n-1\\}$: Locations\n",
    "\n",
    "##### Parameters\n",
    "- $f_{ik}$: Flow between facility $i$ and facility $k$ (represented as `flow[i][k]`)\n",
    "- $d_{jl}$: Distance between location $j$ and location $l$ (represented as `dist[j][l]`)\n",
    "- $n$: Number of facilities (and locations)\n",
    "\n",
    "##### Decision Variables\n",
    "- $x_{ij} \\in \\{0, 1\\}$: Binary variable equal to 1 if facility $i$ is assigned to location $j$, 0 otherwise\n",
    "\n",
    "##### Objective Function\n",
    "\n",
    "$$\n",
    "\\min \\sum_{i=0}^{n-1} \\sum_{j=0}^{n-1} \\sum_{k=0}^{n-1} \\sum_{l=0}^{n-1} c_{ik} \\cdot f_{ik} \\cdot d_{jl} \\cdot x_{ij} \\cdot x_{kl}\n",
    "$$\n",
    "\n",
    "where the coefficient $c_{ik}$ is defined as:\n",
    "\n",
    "$$\n",
    "c_{ik} = \\begin{cases} \n",
    "-1 & \\text{if } (i+k) \\bmod 2 = 0 \\\\\n",
    "1 & \\text{if } (i+k) \\bmod 2 = 1\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "##### Constraints\n",
    "\n",
    "**Assignment constraints - Each facility to exactly one location:**\n",
    "$$\n",
    "\\sum_{j=0}^{n-1} x_{ij} = 1 \\quad \\forall i \\in \\{0, 1, \\ldots, n-1\\}\n",
    "$$\n",
    "\n",
    "**Assignment constraints - Each location receives exactly one facility:**\n",
    "$$\n",
    "\\sum_{i=0}^{n-1} x_{ij} = 1 \\quad \\forall j \\in \\{0, 1, \\ldots, n-1\\}\n",
    "$$\n",
    "\n",
    "**Binary constraints:**\n",
    "$$\n",
    "x_{ij} \\in \\{0, 1\\} \\quad \\forall i, j \\in \\{0, 1, \\ldots, n-1\\}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1565324a",
   "metadata": {},
   "source": [
    "Use this data and the model below to solve the problem. You may notice it take a little while. \n",
    "- If you are running this on your laptop, set $n = 10$ and it should take 20-30 seconds.\n",
    "- If you are running in Colab, set $n = 9$. If it's set to 10, it could take around 3 or 4 minutes. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d911c6c",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 10\n",
    "\n",
    "# Distance between locations\n",
    "dist = [[abs(j - k) + 1 for k in range(n)] for j in range(n)]\n",
    "\n",
    "# Flow/weight between facilities\n",
    "flow = [[(i + k + 1) * ((i - k) ** 2 + 1) for k in range(n)] for i in range(n)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9fcccfd2",
   "metadata": {},
   "source": [
    "Here is the model in `gurobipy`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a41192b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Quadratic assignment problem\n",
    "m1 = gp.Model()\n",
    "m1.Params.Seed = 42\n",
    "\n",
    "# x[i,j] = 1 if facility i is assigned to location j\n",
    "x = m1.addVars(n, n, vtype=GRB.BINARY, name=\"x\")\n",
    "\n",
    "# Build objective function\n",
    "obj = gp.QuadExpr()\n",
    "for i in range(n):\n",
    "    for j in range(n):\n",
    "        for k in range(n):\n",
    "            for l in range(n):\n",
    "                coeff = flow[i][k] * dist[j][l]\n",
    "                # Alternate signs\n",
    "                if (i + k) % 2 == 0:\n",
    "                    coeff = -coeff\n",
    "                obj.addTerms(coeff, x[i, j], x[k, l])\n",
    "\n",
    "m1.setObjective(obj, GRB.MINIMIZE)\n",
    "\n",
    "# Each facility is assigned to exactly one location\n",
    "for i in range(n):\n",
    "    m1.addConstr(x.sum(i, \"*\") == 1, name=f\"assign_facility[{i}]\")\n",
    "\n",
    "# Each location is assigned exactly one facility\n",
    "for j in range(n):\n",
    "    m1.addConstr(x.sum(\"*\", j) == 1, name=f\"assign_to_location[{j}]\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1cc72e6b",
   "metadata": {},
   "source": [
    "Run this without editing any parameters to see how long it takes. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e92b8166",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set parameters here after running in once. \n",
    "\n",
    "# m1.Params.ParameterName = value\n",
    "m1.optimize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "647f576b",
   "metadata": {},
   "source": [
    "Let's try to speed this up. You can...\n",
    "- Try to recall a parameter mentioned in the log file session that may fit this type of problem and searching the [documentation](https://docs.gurobi.com/projects/optimizer/en/current/). This is a good idea. \n",
    "\n",
    "OR\n",
    "\n",
    "- Ask [Gurobot](https://portal.gurobi.com/iam/chat/). Even try giving it the log file just produced to see what it suggests. This is also a good idea. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3079a6d4",
   "metadata": {},
   "source": [
    "<div style=\"\n",
    "  border-left: 6px solid #f59e0b;  /* amber/yellow */\n",
    "  padding:14px 16px;\n",
    "  border-radius:10px;\n",
    "  margin:12px 0;\n",
    "  font-size:16px; line-height:1.4;\">\n",
    "  <strong style=\"font-size:18px;\">🧠 Take a minute to think about this one!</strong><br>\n",
    "  <em>Pause here.</em> There are spoilers just below. \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9dcfeb53",
   "metadata": {},
   "source": [
    "### Click to expand and see the answer"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c52d645",
   "metadata": {},
   "source": [
    "What we are looking for is `PreQLinearize`. Ask [Gurobot](https://portal.gurobi.com/iam/chat/) for an explanation on what this does. To see significant improvement, set `PreQLinearize = 2` in the above a commented out line of code and see what happens. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67617f41",
   "metadata": {},
   "source": [
    "## Handling Infeasibility with IIS and feasRelax"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6851d332",
   "metadata": {},
   "source": [
    "Of course we go back to the good ol' widgets transportation problem. Here is a quick reminder of this now classic (infamous?) problem.\n",
    "\n",
    "We make widgets. Have a set of production facilities that produce boxes of widgets. There is also a set of distribution locations that will then distribute the widgets for sale. Each distribution center has a forecasted demand and each production facility has a min and max number of widgets it can make during this period.\n",
    "\n",
    "But times have changed since we first optimized our original transportation problem. Management has some more strict requirements for production while meeting demand. Each facility must produce at least 90% of its max production. AND the min amount of widgets to be shipped from any production to distribution facility (if it's nonzero) must be at least 50. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa23c0b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "path = 'https://raw.githubusercontent.com/Gurobi/modeling-examples/master/optimization101/Modeling_Session_1/'\n",
    "transp_cost = pd.read_csv(path + 'cost.csv')\n",
    "# get production and distribution locations from data frame\n",
    "production = list(transp_cost['production'].unique())\n",
    "distribution = list(transp_cost['distribution'].unique())\n",
    "transp_cost = transp_cost.set_index(['production','distribution']).squeeze()\n",
    "\n",
    "max_prod = pd.Series([180,200,140,80,180], index = production, name = \"max_production\")\n",
    "n_demand = pd.Series([89,95,121,101,116,181], index = distribution, name = \"demand\")\n",
    "\n",
    "# the min production is a fraction of the max\n",
    "frac = 0.90\n",
    "# min shipment size, if it's > 0\n",
    "C = 50"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f959fedc",
   "metadata": {},
   "source": [
    "### The Original Model\n",
    "\n",
    "- Our **decision variables** are the amount produced at facility $p$ and shipped to distribution center $d$, denoted $x_{p,d}$\n",
    "- We have **constraints** to ensure:\n",
    "  - Min and max production\n",
    "  - Demand is met\n",
    "- The **objective** is to meet demand at **minimal cost**.\n",
    "\n",
    "Our production and distribution `sets` are:\n",
    "$$\n",
    "\\begin{align*}\n",
    "P &= \\{\\texttt{'Baltimore', 'Cleveland', 'Little Rock', 'Birmingham', 'Charleston'}\\}, &\\texttt{production} \\\\\n",
    "D &= \\{\\texttt{'Columbia', 'Indianapolis', 'Lexington', 'Nashville', 'Richmond', 'St. Louis'}\\}  &\\texttt{distribution}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1512496",
   "metadata": {},
   "source": [
    "Model `parameters`:\n",
    "$$\n",
    "\\begin{align*}\n",
    "c_{p,d} &: \\text{cost to ship a widget from} \\space p \\space \\text{to} \\space d &\\texttt{trasp}\\_\\texttt{cost[p,d]}\\\\\n",
    "m_p &: \\text{maximum a production facility} \\space p \\space\\text{can produce} &\\texttt{max}\\_\\texttt{prod[p]}\\\\\n",
    "n_d &: \\text{demand at distribution hub} \\space d &\\texttt{n}\\_\\texttt{demand[d]}\n",
    "\\end{align*}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca626587",
   "metadata": {},
   "source": [
    "As a reminder, here is the original *formulation*:\n",
    "$$\n",
    "\\begin{align*}\n",
    "{\\rm min} &\\sum_{p,d}c_{p,d}x_{p,d}\\\\\n",
    "{\\rm s.t.}\\\\\n",
    "&\\sum_{d}x_{p,d} \\le m_p, &\\forall p \\in P \\quad &\\texttt{can}\\_\\texttt{produce[p]}\\\\\n",
    "&\\sum_{d}x_{p,d} \\ge a*m_p,&\\forall p \\in P \\quad &\\texttt{must}\\_\\texttt{produce[p]}\\\\\n",
    "&\\sum_{p}x_{p,d} \\ge n_d, &\\forall d \\in D \\quad &\\texttt{meet}\\_\\texttt{demand[d]}\\\\\n",
    "&x_{p,d} \\ge 0,  &\\forall p \\in P, d \\in D\\quad &\\texttt{non-negativity}\\\\\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b5e5980",
   "metadata": {},
   "source": [
    "Here is our original model with the additional requirement $x$ is a `semicontinuous` decision variable. Run this scenario as described. \n",
    "\n",
    "I wonder what happens..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6bdfbb2f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# gurobipy code for this above formulation\n",
    "m = gp.Model('widgets')\n",
    "\n",
    "# decision vars\n",
    "x = m.addVars(production, distribution, vtype=GRB.SEMICONT, lb = C, name = 'prod_ship')\n",
    "\n",
    "# constraints\n",
    "can_produce = m.addConstrs((gp.quicksum(x[p,d] for d in distribution) <= max_prod[p] for p in production), name = 'can_produce')\n",
    "must_produce = m.addConstrs((gp.quicksum(x[p,d] for d in distribution) >= frac*max_prod[p] for p in production), name = 'must_produce')\n",
    "meet_demand = m.addConstrs((x.sum('*', d) == n_demand[d] for d in distribution), name = \"meet_demand\")\n",
    "\n",
    "#objective\n",
    "m.setObjective(gp.quicksum(transp_cost[p,d]*x[p,d] for p in production for d in distribution), GRB.MINIMIZE)\n",
    "m.optimize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dfdb62bb",
   "metadata": {},
   "source": [
    "<div style=\"\n",
    "  border-left: 6px solid #f59e0b;  /* amber/yellow */\n",
    "  padding:14px 16px;\n",
    "  border-radius:10px;\n",
    "  margin:12px 0;\n",
    "  font-size:16px; line-height:1.4;\">\n",
    "  <strong style=\"font-size:18px;\">💬 Explainability Moment</strong><br>\n",
    "  <em>This model is infeasible!</em> How would you explain what infeasibility means to a <b>non-technical stakeholder</b> in this context? \n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff3d6ea4",
   "metadata": {},
   "source": [
    "### Computing IIS\n",
    "\n",
    "If you don't know what an IIS is, try and ask [Gurobot](https://portal.gurobi.com/iam/chat)? \n",
    "\n",
    "If you'd rather not do that right now, here's a quick definition.\n",
    "\n",
    "An **IIS (Irreducible Infeasible Subsystem)** is a minimal subset of constraints and variable bounds in an infeasible optimization model that satisfies two key properties:\n",
    "\n",
    "- It is itself infeasible - the subset alone makes the model infeasible.\n",
    "- It becomes feasible if any single constraint or bound is removed - removing any one element from the subset results in a feasible subsystem.\n",
    "\n",
    "And here is how to compute one:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd7558c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "if m.Status == GRB.INFEASIBLE:\n",
    "    m.computeIIS()\n",
    "    \n",
    "    # Check which constraints are in the IIS\n",
    "    for constr in m.getConstrs():\n",
    "        if constr.IISConstr:\n",
    "            print(f\"Conflicting: {constr.ConstrName}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac9cc4a7",
   "metadata": {},
   "source": [
    "<div style=\"\n",
    "  border-left: 6px solid #f59e0b;  /* amber/yellow */\n",
    "  padding:14px 16px;\n",
    "  border-radius:10px;\n",
    "  margin:12px 0;\n",
    "  font-size:16px; line-height:1.4;\">\n",
    "  <strong style=\"font-size:18px;\">💬 Explainability Moment</strong><br>\n",
    "  Now that we see where the issue is, how would you go about fixing this? There are a bunch of ways to do it, so there's no single correct answer. How would you explain each to stakeholders?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a51e688f",
   "metadata": {},
   "source": [
    "Here is the data and model again. Make it feasible! "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6875e982",
   "metadata": {},
   "outputs": [],
   "source": [
    "max_prod = pd.Series([180,200,140,80,180], index = production, name = \"max_production\")\n",
    "n_demand = pd.Series([89,95,121,101,116,181], index = distribution, name = \"demand\")\n",
    "\n",
    "# the min production is a fraction of the max\n",
    "frac = 0.90\n",
    "# min shipment size, if it's > 0\n",
    "C = 50\n",
    "\n",
    "# gurobipy code for this above formulation\n",
    "m = gp.Model('widgets')\n",
    "\n",
    "# decision vars\n",
    "x = m.addVars(production, distribution, vtype=GRB.SEMICONT, lb = C, name = 'prod_ship')\n",
    "\n",
    "# constraints\n",
    "can_produce = m.addConstrs((gp.quicksum(x[p,d] for d in distribution) <= max_prod[p] for p in production), name = 'can_produce')\n",
    "must_produce = m.addConstrs((gp.quicksum(x[p,d] for d in distribution) >= frac*max_prod[p] for p in production), name = 'must_produce')\n",
    "meet_demand = m.addConstrs((x.sum('*', d) == n_demand[d] for d in distribution), name = \"meet_demand\")\n",
    "\n",
    "#objective\n",
    "m.setObjective(gp.quicksum(transp_cost[p,d]*x[p,d] for p in production for d in distribution), GRB.MINIMIZE)\n",
    "m.optimize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9074032",
   "metadata": {},
   "source": [
    "### Using feasRelax\n",
    "IIS computes the subset of constraints in which we know where the problem lies. Another way would be to use `feasRelax`, which relaxes certain parts of the model to make sure it's feasible. In this case we are using `feasRelaxS` which is a little easier to use right out of the box as it relaxes all constraints equally. \n",
    "\n",
    "Below is some code to run a relaxed model, along with what each argument does in the function used. Check out the documentation on [feasRelaxS](https://docs.gurobi.com/projects/optimizer/en/current/reference/python/model.html#Model.feasRelaxS). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec8a9f00",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Relax ALL constraints equally (penalty = 1.0 for each)\n",
    "m.feasRelaxS(0, False, False, True)\n",
    "#            │    │      │      └─ crelax: True = relax constraints\n",
    "#            │    │      └─ vrelax: False = don't relax variable bounds\n",
    "#            │    └─ minrelax: False = just minimize violations\n",
    "#            └─ relaxobjtype: 0=sum, 1=sum of squares, 2=count\n",
    "\n",
    "m.optimize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "713b4177",
   "metadata": {},
   "source": [
    "The model is now minimizing the sum of the total violations. To help us with that, let's first write an `lp file` to see what exactly is happening. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d37a29ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "m.write(\"relax.lp\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed31e4db",
   "metadata": {},
   "source": [
    "The relaxed model adds on a positive and negative violation artificial variable for each constraint. Looking at the lp file, how can we select the artificial variables?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d4abd9a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Analyze violations\n",
    "print(\"\\nViolated constraints:\")\n",
    "for var in m.getVars():\n",
    "    if var.VarName.startswith(  ????  ) and var.X > 1e-6:\n",
    "        print(f\"  {var.VarName}: {var.X:.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24418a6f",
   "metadata": {},
   "source": [
    "So there is a problem with the production amount in a particular production location. Specifically, the minimum amount we are *requiring* to be produced. Go back to the lp file to investigate further. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aaebab37",
   "metadata": {},
   "source": [
    "Let $A_{\\hat{p}}$ be the artificial variable noted above for production location $\\hat{p}$. The constraint(s) in question, looking at the lp file:\n",
    "\n",
    "$$\n",
    "\\begin{align*} \n",
    "\\sum_{d}x_{\\hat{p},d} + A_{\\hat{p}} &\\ge a*m_p\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "912cd7e8",
   "metadata": {},
   "source": [
    "<div style=\"\n",
    "  border-left: 6px solid #f59e0b;  /* amber/yellow */\n",
    "  padding:14px 16px;\n",
    "  border-radius:10px;\n",
    "  margin:12px 0;\n",
    "  font-size:16px; line-height:1.4;\">\n",
    "  <strong style=\"font-size:18px;\">💬 Explainability Moment</strong><br>\n",
    "Now that you see what this looks like algebraically, how would you communicate this?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "191cbdda",
   "metadata": {},
   "source": [
    "<div style=\"\n",
    "  border-left: 6px solid #3b82f6;  /* blue */\n",
    "  padding:14px 16px;\n",
    "  border-radius:10px;\n",
    "  margin:12px 0;\n",
    "  font-size:16px; line-height:1.4;\">\n",
    "  <strong style=\"font-size:18px;\">📘 Optional Activity: Try different argument values for feasRealaxS</strong><br>\n",
    "  For example, we minimized the sum of total constraint violations. Check the other options and see how that changes things. \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70540c9d",
   "metadata": {},
   "source": [
    "#### Homework problem!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1e06249",
   "metadata": {},
   "source": [
    "<div style=\"\n",
    "  border-left: 6px solid #3b82f6;  /* blue */\n",
    "  padding:14px 16px;\n",
    "  border-radius:10px;\n",
    "  margin:12px 0;\n",
    "  font-size:16px; line-height:1.4;\">\n",
    "  <strong style=\"font-size:18px;\">📘 Optional Activity: Use your new skills!</strong><br>\n",
    "  The below model is infeasible. Fix it however you see fit!\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0236520b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# gurobipy code for this above formulation\n",
    "m2 = gp.Model('widgets')\n",
    "\n",
    "max_prod2 = pd.Series([210,225,140,130,220], index = production, name = \"max_production\")\n",
    "n_demand = pd.Series([89,95,121,101,116,181], index = distribution, name = \"demand\")\n",
    "\n",
    "frac = 0.75\n",
    "# decision vars\n",
    "x2 = m2.addVars(production, distribution, vtype=GRB.SEMICONT, lb = C, name = 'prod_ship')\n",
    "y2 = m2.addVars(production, vtype=GRB.BINARY, name = 'prod_on')\n",
    "\n",
    "# constraints\n",
    "can_produce = m2.addConstrs((gp.quicksum(x2[p,d] for d in distribution) <= max_prod2[p] for p in production), name = 'can_produce')\n",
    "must_produce = m2.addConstrs((gp.quicksum(x2[p,d] for d in distribution) >= frac*max_prod2[p] for p in production), name = 'must_produce')\n",
    "meet_demand = m2.addConstrs((x2.sum('*', d) == n_demand[d] for d in distribution), name = \"meet_demand\")\n",
    "xy_link = m2.addConstrs((x2[p,d] <= max_prod2[p]*y2[p] for p in production for d in distribution), name = 'xy_link')\n",
    "only_four = m2.addConstr((y2.sum() <= 4), name = 'only4')\n",
    "\n",
    "#objective\n",
    "m2.setObjective(gp.quicksum(transp_cost[p,d]*x2[p,d] for p in production for d in distribution), GRB.MINIMIZE)\n",
    "m2.optimize()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "gurobi_ml",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
